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The possibility of proton acceleration to very high energies in astrophysical sources may have 
unique observational consequences. In particular, the decay of secondary mesons created in the 
interaction of energetic protons with photons or nucleons gives rise to high-energy gamma rays 
and neutrinos with potentially observable fluxes. Recently, there has been considerable interest in 
the parameterization of the energy spectra of these secondaries. Less attention has been paid to 
the angular distributions, which may have an important effect on observational quantities and are 
required to address collisions between protons with different energies and an arbitrary scattering 
angle. In this work, we study the complete particle distributions of secondary mesons created in 
proton - proton collisions. We present parameterizations of the energy and rapidity distributions 
of secondary pions and kaons that reproduce results generated with the event generator PYTHIA 
to within ~10% in the bulk of the parameter space. The parameterizations are based on incident 
proton energies from 1 TeV to 1 PeV and are suited for extrapolation to higher energies. We 
present several applications of these parameterizations. Energy spectra and angular distributions 
of stable decay products (electrons, neutrinos and gamma rays) follow readily. We give an example 
of the gamma-ray spectrum that results from the decay of 7r° mesons created in a proton - proton 
collision. We show that there is a strong correlation between the energy of secondary mesons and the 
degree of collimation around the direction of the colliding protons. This effect may have important 
implications for the detection possibility of neutrinos created in the interaction of a developing GRB 
with its surroundings. 

PACS numbers: 13.60.Le, 25.40.Ep, 13.85.-t, 98.70.Sa 



I. INTRODUCTION 

The possibility of proton acceleration to very high 
energies in astrophysical sources may provide unique 
observational opportunities. The interaction of ener- 
getic protons with photons or nucleons results in copi- 
ous production of secondary mesons decaying into high- 
energy gamma rays and neutrinos that can be observed 
with current and future detectors. The recently ob- 
served TeV gamma-ray emission from supernova remnant 
RX J1713.7-3946 [l| has been attributed to this mecha- 
nism [2] , although such an origin is still under debate Q . 
TeV gamma rays have been reported in coincidence with 
gamma-ray burst (GRB) 970417a 0, but also in this 
case it is not established whether the origin is hadronic 
i- 

The existence of astrophysical proton accelerators is 
indicated by observations of high-energy cosmic rays 
(CRs). There is evidence for a substantial proton compo- 
nent above the 'knee' at ~4 x 10^ GeV in the cosmic-ray 
spectrum (see e.g. Ref. Q for a review). Observations 
of extensive air showers due to CRs with energies up 
to ~10^^ GeV are consistent with nucleon primaries, al- 
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though other primaries are also possible (e.g., Ref. 8]). 

Various astrophysical systems have been suggested as 
CR sources. Galactic supernova remnants are the lead- 
ing candidate for the generation of GRs with energies up 
to ~10* GeV (see e.g. Ref. ^Qj). Several extragalactic 
sources have been considered as possible sources of higher 
energy GRs, such as active galactic nuclei (see how- 
ever Ref. Ifljl), hot spots of Fanaroff-Riley class II radio 
galaxies [unil, pulsars 13] and GRBs 

A population of high-energy protons in these sources 
would carry a rich phenomenology. In GRBs for example, 
the interaction of accelerated pro tons with GRB photons 
leads to - 10^ GeV neutrinos ^ and to -10^ - 10^ GeV 
gamma rays [l3] • High-energy proton interactions may 
play an important role in the interaction of a developing 
GRB with its environment, e.g. when the fireball has 
not yet emerged from the stellar surface [H, [ll| or when 
energetic GRB protons collide with cold protons in the 
GRB surroundings 0,[2l|,[22. 

Detailed parameterizations of the energy spectra of 
secondary particles created in proton - proton (pp) col- 
lisions are essential in the study of particle production 
in astrophysical proton accelerators. Such parameteriza- 
tions were recently presented by Kamae et al. [l^ and 
Kelner, Aharonian and Bugayov 24]. However, the pa- 
rameterizations presented by these authors do not in- 
clude the angular distributions of the secondary parti- 
cles. As a consequence, these parameterizations can only 
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be applied to the scattering geometry for which they were 
derived, viz. with a target proton at rest. This is suffi- 
cient if one of the protons is at rest in the observer frame, 
but a more general treatment is required if the bulk of 
the protons can be accelerated. 

The parameterization of the complete particle distri- 
butions is an important generalization because this pro- 
vides, through Lorentz transformations, secondary par- 
ticle distributions and energy spectra after a collision of 
two protons with arbitrary energies and an arbitrary in- 
cident angle. Such a parameterization can therefore be 
applied to any acceleration scenario. A second advantage 
is that the full distribution includes correlations between 
the energy and the angle of outgoing particles, which 
may have important effects on observable quantities in a 
non-isotropic environment. 

Badhwar, Stephens and Golden [2^ , Stephens and 
Badhwar ^201 and Blattnig et al. [27'| have presented pa- 
rameterizations of the complete distributions of charged 
and neutral pions and charged kaons created in pp col- 
lisions. However, these parameterizations are valid for 
incident proton energies Ep < 2 x 10^ GeV, which is 
much lower than the highest proton energies ~10^^ GeV 
expected in accelerating astrophysical sources. 

In this paper we study the complete distributions of 
secondary particles produced in pp collisions through 
Monte Carlo simulations. We consider a proton with en- 
ergy 10^ GeV < Ep < 10^ GeV that coUidcs with a pro- 
ton at rest, which corresponds to center-of-mass energy 
43 GeV < < 1.4 X 10^ GeV. We assume that the dis- 
tribution of a secondary particle species is invariant under 
rotations around the collision axis, which implies that the 
distribution is fully parameterized with two independent 
kinematical variables. We present parameterizations of 
the energy and rapidity distributions of secondary pions 
and kaons. The parameterizations are based on Monte 
Carlo data in the simulated energy range but can be ap- 
plied to collisions of protons with higher energies. 

We consider only secondary pions and kaons and not 
their stable decay products, viz. electrons, neutrinos and 
gamma rays. This approach separates the physics in the 
pp collision from subsequent decay processes. Energy 
spectra and particle distributions of the resulting sta- 
ble daughter particles are readily found from our results 
(either analytically or as part of a computer code) and 
the well-known decay spectra of pions and kaons (see e.g. 
Refs. ^,231)- We do not separately consider short-lived 
mesons (such as 77, p or uj) because their lifetime is much 
shorter than that of charged kaons and pions. The decay 
products of these mesons, mostly pions and gamma rays, 
are grouped together with the prompt secondaries. The 
restriction to the pp interaction per se gives our results a 
broad applicability. For example, it has been pointed out 
recently [2^[^ that energy losses of pions and kaons can 
leave an imprint on the energy spectra of the daughter 
particles in GRB jets. A proper treatment of this effect 
requires knowledge of the pion and kaon distributions. 

Hadron interactions have a complex phenomenology 



due to the compositeness of the ingoing and outgoing par- 
ticles. It is currently not possible to compute the cross 
section or the resulting particle distribution from first 
principles. Therefore, detailed studies of particle produc- 
tion in hadron interactions require the use of event gen- 
erators. In this work, we simulate the pp interaction with 
the event generator PYTHIA [sTj, which is tested against 
experimental data. It is capable of simulating various 
incident and target particles so that it is possible to ex- 
tend this work to proton - neutron and proton ~ photon 
interactions with essentially the same code. PYTHIA uses 
the 'Lund string model' 32] to describe the process of 
secondary meson formation. 

This paper is organized as follows: in section |TT1 we 
present experimental data on the pp cross section and the 
charged multiplicity, i.e. the number of charged particles 
created in a single inelastic collision. In section IIIII we 
discuss the kinematics of the simulated interaction and 
introduce the particle distribution with respect to en- 
ergy and rapidity. Details on the event simulation with 
PYTHIA and the fitting procedure are discussed in sec- 
tion |TVl In section |Vl we present a comparison between 
PYTHIA results and experimental data and we present the 
parameterizations of the energy spectra and particle dis- 
tributions of secondary pions and kaons. Applications 
of these parameterizations are considered in section IVII 
We demonstrate through explicit examples how the pa- 
rameterizations can be used to study particle production 
in collisions of protons with different energies and an ar- 
bitrary incident angle. We also present an example in 
which we derive the gamma-ray energy spectrum result- 
ing from tt'^ 77 decay. In section IVIIl we discuss the 
application of the parameterizations to incident protons 
with very high energies. We discuss the results in section 
IVIIII Conclusions are presented in section HXl 



II. EXPERIMENTAL DATA ON THE CROSS 
SECTION AND SECONDARY MULTIPLICITY 
IN PROTON - PROTON INTERACTIONS 

In this section, we discuss experimental data on the 
cross section of pp interactions and on the number of 
charged particles created in a pp interaction. The data 
presented in this section are used in section |V] to validate 
our numerical method. 



A. Cross section 

Proton - proton interactions are usually separated 
into elastic scattering, in which no particles are created; 
diffractive interactions, in which the energy transfer be- 
tween the protons is small; and inelastic non-diffractive 
interactions, in which the energy transfer is large enough 
for the constituent quarks and gluons to interact. The 
total pp cross section dtot can be expanded in terms of 
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FIG. 1; Comparison between the inelastic pp cross section 
calculated with PYTHIA (open squares) and experimental data 
(disks). The solid line represents the fit given in eq. ((3}. 
Experimental data is taken from the PPDS (see footnote [T|. 
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FIG. 2: Comparison between the charged multiplicity cal- 
culated with PYTHIA (squares) and experimental data (other 
symbols). Open symbols correspond to inelastic processes 
(inel), solid symbols correspond to non-single-diffractive 
(NSD) processes. The solid line shows the approximation 
given in eq. ((5)). 



section crinei = <J-ad + Csd + Cdd becausB all processes that 
create secondary particles are contained in this quantity. 
At energies above the threshold energy £^th = 1-22 GeV 
and below ^/s = 3 x 10^ GeV, the inelastic cross section 
for a proton with energy Ep interacting with a target 
proton at rest can be fitted with (U '■ 

'^fLi[Ep) - (33.24 - 3.624 log .Ep + L325 [logEpf 

X I 1 - ( ^ I 1 mb , (2) 

where Ep is measured in GeV. In deriving this for- 
mula, it is assumed that the ratio of the inelastic cross 
section to the total cross section, which for energies 
43 GeV < Vs < 63 GeV is given by [H 




CTtot = T21o- 



inel : 



(3) 



holds for higher energies 63 GeV < ^/s < 3 x 10^ GeV as 
well. The incident proton energy Ep in eq. ([2]) is related 
to the center-of-mass energy ^/s as 



E.. 



^ 2mpC^ 



(4) 



where nip is the proton mass. In figure [1] we show the 
approximation given in eq. ([2]) together with PYTHIA re- 
sults (see section|V]below) and the available experimental 
data.^ This shows the validity of approximation 



B. Secondary multiplicity 

Bubble chamber and accelerator experiments have 
shown that the number of charged particles created in 
proton - (anti)proton collisions, i.e. the charged mul- 
tiplicity, increases as a function of the incident proton 
energy.^ We find that up to the highest energies cur- 
rently accessible, ^/s < 1.8 x 10^ GeV, the charged par- 
ticle multiplicity in non-diffractive pp interactions is well 
fitted with 



M«i5(s) =0.89 



1.24 logs 



0.34 log^s + 0.077 log^s 



,(5) 



these processes as 

O'tot = CTnd + CTsd + (Tdd + CTcl , (1) 

where Und , Csd , crdd and del are the cross section for non- 
diffractive processes (the hard QCD processes), single 
diffraction {AB -> XB or AB -> AX), double diffrac- 
tion {AB — > XY), and elastic scattering {AB AB) 
respectively. In this work we do not explicitly separate 
diffractive and non-diffractive processes because we are 
mostly interested in astrophysical applications where it 
will be impossible to distinguish between these compo- 
nents; see Ref. [2^ for a separate treatment. 

We are primarily interested in the inelastic pp cross 



This functional form is an extention of an approximation 
due to Matthiae [36| (see also Ref. [13]) which is valid 
only up to y/s < 540 GeV. The last logarithmic term, 
which does not appear in the approximation by Matthiae, 



^ A compilation of experimental data on the inelastic pp cross sec- 
tion is available at the Particle Physics Data System (PPDS) 
website http: //wwwppds . ihep . su: 8001/. In producing figure [Tl 
we have not considered experimental data that is marked with 
the warning comment 'W. 

^ The charged particle multiplicity in pp and pp interactions is 
virtually identical at ISR energies = 53 GeV (Ref. [34l |: see 
also Ref. [H). 
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is required in order to fit the multiplicity at both high 
and low energies. We present in figure [5] experimental 
data"^ together with the approximation given in eq. ([5|) 
and PYTHIA results (see section IVl below) . 

A logarithmic dependence A^ch oc log s is commonly 
interpreted as due to an increase in phase space because 
the range of allowed rapidities scales as log(s/m^c^) (e.g., 
Ref. ^35.]). A stronger increase in secondary multiplicity 
is then attributed to an additional rise in the level of the 
observed central rapidity plateau, the origin of which is 
not understood from first principles. 

At high energies, data on the neutral particle multi- 
plicity is sparse because of experimental difficulties. As 
a result, there is no fit to the neutral particle multiplicity 
that extends to ^/s > 50 GeV which is based on experi- 
mental data. A fit to the separate multiplicities of both 
charged and neutral pions and charged kaons created 
in pp collisions for center-of-mass energies < 53 GeV 
was presented in Ref. [4^ . 

The scarcity of experimental data on separate parti- 
cle multiplicities at high energies motivates the use of 
event generators such as PYTHIA. In section IV Al we 
show that PYTHIA correctly reproduces experimental re- 
sults on the total charged multiplicity. In section IV Bl 
we present a fit to PYTHIA results on charged and neu- 
tral pion and kaon multiplicities in the energy range 
43 GeV < < 1.4 X 10^ GeV. 



For given particle energy e, the rapidity cannot take any 
value. The mass-shell relation implies that —yi < y < yi, 
where 



yi = arccosh ( - 



(7) 



and m is the secondary (pion or kaon) mass. A second 
requirement follows from energy conservation in the pp 
collision. If the energy of the secondary particle e > 
iTipC^, the rapidity is additionally bounded hy y > y2, 
where 



2/2 = arctanh ( — 



(8) 



In this equation (3'^ is the proton velocity in the center- 
of-mass frame in units of c, which we take to be equal 
to one for incident proton energies Ep ^ rUpC^ in the 
following calculations. Note that eq. ^ can be applied 
in any frame, while eq. ([8]) only holds in the lab frame. 



B. Secondary particle distribution 

We are interested in the particle distribution for one- 
pion and one-kaon inclusive processes. 



pp 



XY, 



(9) 



III. KINEMATICS AND SECONDARY 
PARTICLE DISTRIBUTION 

In this and the following sections, we consider an en- 
ergetic proton that moves along the z-axis and collides 
with a proton at rest, i.e. a fixed target. This scattering 
geometry is referred to as the lab frame. We use pz to 
denote a longitudinal momentum, along the z-axis, and 
Pt to denote a transverse momentum. 



A. Kinematics 

Assuming that the secondary particle distribution is 
symmetric around the collision axis, the phase space of 
the outgoing particles is fully parameterized with two 
independent kinematical variables. Here, we choose the 
energy e and the rapidity y, which is defined as 



V = t: In f — 



PzC 



PzC 



tanhy = 



PzC 



(6) 



where X denotes a single pion or a single kaon and Y 
may be any combination of particles with the appropriate 
quantum numbers. We denote by n(e, y)dedy the number 
of created particles of a given species with energy and 
rapidity in the range (e . . . e -I- de) x (y . . .y + dy) : 



n{e,y) 



d^N 



1 



dedy aind dedy ' 



(10) 



where CTinci = fnd + Csd + Cdd is the inelastic pp cross sec- 
tion and a is the inclusive cross section to detect a par- 
ticle of a given kind (assuming an ideal detector). This 
cross section is equal to the weighted sum of n-particle 
exclusive cross sections cr„ (i.e., the cross section to create 
exactly n particles)'*: 



a = ^nan = Alcrinol : 



(11) 



where M = M{s) is the multiplicity of the given particle 
species. The particle distribution n(e, y) is related to 
the Lorentz invariant differential cross section ed^a/dp^, 



3 The experimental data is taken from Refs. [37l.l3i.[39l . l40 l l4ll. l4^ . 

The data at ^ = 1.8 X 10^ GeV is obtained by the E735 ex- 
periment [43II . which does not cover the full particle phase space. 
We use results from Ref. |40|| . who have determined the charged 
particle multiplicity through a fit to experimental data on the 
multiplicity distribution. 



* We do not consider the exclusive cross sections separately be- 
cause we are interested in particle creation by all processes to- 
gether. To a first approximation, the relative sizes of the n- 
particle exclusive cross sections depend on energy only through 
the total multiplicity [45l |. This 'KNO scaling' is known to be 
violated at energies > 500 GeV [sllil. 
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which is often used to represent experimental data, as 
follows: 

nie,y)^ — [ ^ U ,, , ■ 12 

IV. NUMERICAL METHOD 

A. Configuration of tiie PYTHIA event generator and 
initial conditions 

The pp interaction is simulated with PYTHIA version 
6.324 using default values for most of the control param- 
eters. Elastic and diffractive processes are included by 
selecting MSEL=2. In comparing PYTHIA results to ex- 
perimental data on the cross section and charged multi- 
plicity, we allow for pion and kaon decay. In determin- 
ing the parameterizations of the particle distributions, 
pion and kaon decay are switched off with the command 
MDCY(PYCOMP(ID) ,1)=0, where ID is the corresponding 
particle identification number. This approach separates 
the physics in the pp collision from subsequent processes, 
such as secondary synchrotron emission prior to decay, 
etc. 

The PYTHIA code relies on the Lund string model [s^l 
to determine the distribution and multiplicity of sec- 
ondary particles. We present a brief discussion of this 
model in appendix |^ more details can be found in the 
PYTHIA manual [sij . 

We simulate pp collisions for incident proton energies 
Ep = 10^ GeV, Ep = 10^ GeV, Ep = 10^ GeV and 
Ep — 10^ GeV colliding with a proton at rest. For higher 
values of the incident proton energy, PYTHIA signals a 
loss of accuracy in kinematical variables in some of the 
generated events. 

B. Fitting procedure 

The secondary particle distributions are discretized, 
spanning the full range of available energy and kinemat- 
ically allowed rapidity. In this process, the energy is 
divided into 200 bins with size Aci with a logarithmic 
division and the rapidity is divided into 100 bins with 
size Ajji with a linear division. The logarithmic energy 
division is chosen because we consider up to seven energy 
decades; the rapidity division is linear because the range 
of allowed rapidities scales with the logarithm of energy. 
The number of bins is limited by computational issues, as 
data files become increasingly large and fitting becomes 
increasingly time-consuming with an increasing number 
of bins. We have verified that this number of bins is suf- 
ficient for convergence of the resulting parameterization. 

We use MINUIT^ as a minimization algorithm for the 



^ CERN Program Library entry D506; doc- 



weighted squared difference between the PYTHIA results 
and the particle distribution fit function n(e,y).^ We 
consider only statistical errors in the PYTHIA results. We 
simulate iVov = 10^ collisions for every incident proton 
energy, which results in a statistical error of a few per- 
cent near the maximum values of en{e,y). We compare 
our results with parameterizations based on other event 
generators to obtain an estimate of the importance of 
systematic uncertainties within the models underlying 
PYTHIA in section IVlira 

The relative deviation between a PYTHIA data point Ui 
and the fitted value n(ei, yi) is expressed as 

S, = ^ ' , (13) 

n^ 

where rii is the number of particles in a bin with average 
energy and average rapidity yi divided by the bin size 
Aei X Ayi. We note that the deviations are expected to 
follow a Gaussian distribution with average value 

(6.) cx (14) 

where the dependence on energy is due to the logarithmic 
energy binning. In particular, the average deviation size 
is expected to be roughly independent of energy for a 
n(e) cx energy spectrum. 



V. RESULTS 

A. Comparison of PYTHIA results with experimental 

data 

We show in figures [T] and O experimental data on the 
pp cross section and charged multiplicity together with 
PYTHIA results. In producing these figures, we have not 
switched off any natural particle decays in the PYTHIA 
simulations in order to compare PYTHIA results with ex- 
perimental data (cf. section HV Ap . 

We observe from figure [1] that the PYTHIA cross sec- 
tions are compatible with an extrapolation of the ex- 
perimental data (see footnote [T]) in the energy range 
43 GeV < < 1.4 X 10^ GeV. In figure El we show a 
comparison between experimental data and PYTHIA re- 
sults on the charged multiplicity. In this work, we are 
interested in particle creation by all inelastic processes. 
However, experimental data on the charged multiplic- 
ity resulting from all inelastic processes is available only 
up to = 62 GeV [H, |4l|, SI]. Experimental data 
on the charged multiplicity resulting from the restricted 



umentation is available on the website 

http: //wwwasdoc .web . cern. ch/wwwasdoc/minuit/ 

We do not explicitly write the dependence of n on Ep here and 
in the following sections. 
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TABLE L Average fraction of the incident proton energy car- 
ried by the outgoing particle species. 

Np Ns 7 7r+ TT" 7r° K+ K' K° 'K^ 
f 0.56 0.033 0.013 0.13 0.095 0.12 0.016 0.011 0.014 0.011 



TABLE IL Numerical values for the constants in the multi- 
plicity approximation formula (|15|l . 





n+ 






K+ 


R- 


K° 




CO 


4.5 


3.8 


4.9 


0.49 


0.32 


0.36 


0.29 


Cl 


-1.7 


-1.7 


-2.1 


-0.23 


-0.20 


-0.17 


-0.18 


C2 


0.50 


0.50 


0.58 


0.063 


0.060 


0.054 


0.054 



class of non-single-diffractive (NSD) interactions is avail- 
able up to much higher energies y/s ~ 1.8 x 10'^ GeV 
[13, Eo, El, 1131 ■ To verify our numerical method, we have 
performed a separate simulation^ of NSD interactions to 
compare the NSD charged multiplicity with experimental 
data. We show in figure [2] that the PYTHIA results on the 
charged multiplicity due to both inelastic processes and 
NSD processes are compatible with experimental data. 

B. Average secondary energy and multiplicity 

The fraction / of the incident proton energy carried by 
a certain secondary particle species is virtuall y in depen- 
dent of the incident proton energy (see Ref. [4^). The 
average fractions for nucleons, photons, pions and kaons 
found by PYTHIA simulations are given in table [J In this 
table, Np and Ng denote primary and secondary nucle- 
ons, respectively (see below). Other possible secondaries 
(direct electrons, muons, neutrinos) together carry less 
than 0.1% of the incident proton energy. 

We define the primary nucleon as the most energetic 
outgoing nucleon. The probability that the primary nu- 
cleon is a proton is 0.70; if it is a proton, it carries an 
average fraction 0.63 of the incident proton energy. The 
probability that the primary nucleon is a neutron is 0.30; 
if this is the case, the average energy fraction is 0.41. The 
energy fraction carried by the primary nucleon as shown 
in table H] represents the weighted average. 

We fit PYTHIA results on the secondary par- 
ticle multiplicities within the energy range 
43 GeV < < 1.4 X 10^ GeV. We find that both 
charged and neutral pion and kaon multiplicities are well 
approximated with the following function: 

Ml = Cq+Ci log S + C2 log^ s , (15) 



For the NSD case, we have switched off single diffraction 
in the event generation with the commands MSUB(92)=0 and 
MSUB(93)=0. 



PYTHIA results 




e (GeV) 



PYTHIA results 
Fit 




10° lO' 10^ 10-' lO'' 10^ lO'' 
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FIG. 3: Energy spectra of vr^ (top panel) and (bottom 
panel) mesons created in a collision with incident proton en- 
ergy Ep = 10® GeV. Note the different scales on the vertical 
axes. Upper graphs: comparison of PYTHIA results and fit to 
the energy spectrum normalized as en(e); lower graphs: devi- 
ation Si — n{ei,yi)/ni — 1 between PYTHIA results points and 
fitted values. 



where cq, ci and C2 are numerical constants whose values 
are given in table |TT] and s is expressed in units of GeV^ . 
The charged kaon multiplicity deduced from eq. (jlSp 
is within ~5% of experimental data at ^/s = 45 GeV 
(43 |. The charged pion multiplicities determined by this 
equation are ^^10% lower than the experimental values. 
This discrepancy can be attributed to the fact that we 
have considered only prompt pions (i.e., excluding pions 
from kaon decay). 
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C. Pion and kaon energy spectra 

We present in figure [3] the 7r+ and energy spectra 
resulting from a collision of a proton with incident en- 
ergy Ep = 10^ GeV with a proton at rest. We find that 
the energy spectra for all secondary particles and inci- 
dent proton energies 10'^ GeV < Ep < 10^ GeV are simi- 
lar in shape. To a first approximation, the energy spec- 
tra follow a power-law, reflecting the absence of an 
energy scale between the secondary mass and the maxi- 
mum available energy. This is supplemented with addi- 
tional functions that we denote with a{e), P{e), 71(e) and 
72 (e) (here and in the following we do not explicitly write 
the dependence of the model parameters on the incident 
proton energy Ep to avoid cluttering of the notation). 
Thus we write the pion and kaon energy spectrum in the 
following way: 

n(e) = noe-ia(e)/3(e)7i(e)72(e) , (16) 

where no is a normalization constant, a(e) accounts for 
the convex shape on a log-log scale, /3(e) incorporates an 
exponential decline at higher and lower energies, 71(e) is 
a strong cutoff near the mass threshold and 72(e) is a 
strong cutoff near the maximum available energy. These 
functions are parameterized as follows: 



no 


= 1.21 X IQPo+p^pI ■ 


(17a) 


a(e) 


— ,Pl(log(c)-2p2) . 


(17b) 


/3(e) 


= 10"^'/'=*^"'° 10"^'/'=^^"™ ; 


(17c) 


71(e) 


= tanh (prolog (e/mc^)) ; 


(17d) 


72(e) 


= tanh(p8olog(i;p/e)) , 


(17e) 



where €3 = 10^"+^*, £5 = IQp^+pb and all energies are 
expressed in units of GeV. The following parameters vary 
with incident proton energy: 

Po Poo +Poilog(£'p) ; (18a) 

Pi = Pw + Pii log(£'p) ; (18b) 

P2 P20 + P21 log(£'p) ; (18c) 

Pi = P40 - P2 ; (18d) 

P6 = Pea +P2 ■ (18e) 

Thus, the energy spectrum of secondary pions and kaons 
is fully described in terms of 12 free parameters pij for 
every particle species. These parameters and their nu- 
merical values, which are determined by a least-squares 
fit, are given in table IIIll 

For pions, deviations between fit values n{ei) and 
PYTHIA results are less than 5% except for very high 
energies (e > Ep/2) and some occasional points near 
the mass threshold where the deviation is ~10%. For 



Figures of those and others fits are available at 
http : //www.nikhef .nl/~hkoers/ppf it. 



kaons, statistical fluctuations are larger since the number 
of kaons to pions is roughly 1:10. At intermediate ener- 
gies the fit is nevertheless within ^--^5% of PYTHIA results 
except for some isolated points. Near the mass thresh- 
old deviations increase to ~20%; at very high energies, 
where en{e,y) is typically more than an order of mag- 
nitude smaller than its maximum value, deviations can 
increase to ^--^40%. We have verified that the parameter- 
ized spectra integrate to the right multiplicities as given 
in eq. (fT5|) within a few percent, except for the K° spec- 
trum for which the deviation is ~10% at the low end of 
the simulated proton energy range. 



D. Pion and kaon energy and rapidity distributions 

We present the pion and kaon rapidity distributions, 
i.e. n(e,y) at fixed e = eo, for incident proton en- 
ergy Ep — 10^ GeV and secondary particle energy 
e = lO'^ GeV in figured! We find that rapidity distribu- 
tions for different proton energies and different secondary 
particle energies are very similar in shape. This shape is 
different for pions and for kaons, hence we treat pions 
and kaons separately in the following. 



1. Pions 

The pion rapidity distributions at fixed energy are 
found to be approximately Gaussian near their maxi- 
mum values (see fig. HJ. At intermediate pion energies, 
e ^ -y/ Ep, the distributions exhibit a low-rapidity tail 
that falls off exponentially (all energies are expressed in 
units of GeV). The distributions fall off very steeply at 
the boundaries of the kinematical domain given in eqs. 
and JH]). 

We factorizc the full particle distribution n(e, y) into a 
modified energy spectrum n(e) and a rapidity-dependent 
function (/)(e, y) that contains both a Gaussian and an 
exponential part: 

n^{e,y) = h{e)(p{e,y) , (19) 

where 

n(e) = noe-ia(e)/3(e)72(e)10^°+293 . (20a) 

(t)ie,y) = lO-V939?(y-92)H9i ^ (20b) 

and no, a(e), /3(e) and 72(e) are defined in eqs. (|17p . The 
parameters qi depend on the pion energy e and on the 
incident proton energy Ep in the following way (here and 
in the following we do not explicitly write the dependence 
of the parameters qi on the pion and proton energies to 
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FIG. 4: Rapidity distributions of 7r+ (top panel) and K'^ 
(bottom panel) secondaries created in a collision of a proton 
with energy Ep = 10^ GeV with a proton at rest. The sec- 
ondary particle energy is taken to be e = 10^ GeV. Upper 
graphs: comparison of PYTHIA results and fit to the rapidity 
distribution en{e,y) as a function of y; lower graphs: devi- 
ation Si = n{ti,yi) /rii — 1 between PYTHIA results and fitted 
values. 



avoid cluttering of notation): 

go = 900 + 901 ; 

qi = 910 + 911 + 912)^ + 913^^ ; 

92 = ln(e) + 920 + 921 (log(e) + 1) 

20''23 log(-Ep) log(m/e) . 

93 = 930 + 10«^^+«^^'°^^'\ 

where we introduced the variable 
C = 21og(e)/log(i?p) 



(21a) 
(21b) 

J^Q922log(Ep)log(e/Bp) 

(21c) 
(21d) 



1 



(22) 



every pion species (tt^, tt~ and 7r°). The fitted values for 
the coefficients are given in table IIVI 

Deviations between the parameterizations and PYTHIA 
results are within 10% in the range in which the rapidity 
distribution is within one order of magnitude of the max- 
imum value and for pion energies 1 GeV < e < 0.1 i?p, 
except for a few isolated points that are typically within 
20%. At high energies, e > 0.1 Ep, deviations increase to 
~30% at the borders of the considered rapidity interval, 
in concordance with eq. . 

We have verified that integrating the energy and ra- 
pidity distributions over rapidity reproduces the energy 
spectra. The deviations between these spectra and 
PYTHIA results are similar to the deviations for the di- 
rect fit to the energy spectra (see section IV C|) . except 
at very low energies e < 2m7r where deviations increase 
to ~30%. The multiplicities obtained by integrating the 
distributions over energy and rapidity are within a few 
percent of those given by eq. ([T5]). 



2. Kaons 

The shape of the kaon rapidity distributions is similar 
to the low-rapidity part of the pion rapidity distributions 
(see fig. U). We find that the kaon energy and rapidity 
distributions are well described with: 



(23) 



where (j}{e, y) is defined in terms of model parameters qi 
in eqs. (|20p and n(e) is a modified energy spectrum: 



n(e) = noe-^a(e)/3(e)10*+2«=' . 



(24) 



Hence, we have parameterized the pion energy and ra- 
pidity distributions in terms of 24 free parameters qij for 



The quantities no, a(e) and /3(e) are defined in eqs. (fT7)) . 
We find that the parameterizations for qi given in eqs. 
(Plj) approximate the PYTHIA results well if we fix 523 = 0. 
Therefore, the kaon energy and rapidity distributions are 
fully parameterized in terms of 23 free parameters for 
every kaon species {K^ , , and K'^). The fitted 
values for these parameters are presented in table IIVI 

Deviations between the approximation (|23p and 
PYTHIA results are similar to the deviations for the pa- 
rameterizations of the pion distributions, except that 
fluctuations are larger. This results in deviations up to 
~3G% at isolated points for all energies. 

For the and mesons, integrating the full distri- 
butions over the rapidity reproduces the energy spectra 
with deviations similar to those for the direct parameter- 
izations of the energy spectra presented in eq. p6p . For 
and mesons, deviations near the mass threshold 
are ~30%, while at very high energies (e > Ep/2) the 
deviations can increase to ^--^50%. These large deviations 
occur only at energies where e n(e, y) is more than an or- 
der of magnitude smaller than the maximum value. The 
multiplicities obtained by integrating the parameterized 
distributions (j23p over energy and rapidity are within a 
few percent of the values given by eq. (E 
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FIG. 5: The tt" distribution en{e,y) as a function of energy 
e and rapidity y after a collision of a proton with energy 10® 
GeV with a proton at rest (lab frame). The discretization 
and the observed 'floor' are for presentational purposes. 



FIG. 6: The tt" distribution en{e,y) as a function of energy 
e and rapidity y after a collision of two protons with energy 
730 GeV (center-of-mass frame). The discretization and the 
observed 'floor' are for for presentational purposes. 



VI. APPLICATIONS 

In this section, we consider applications of the param- 
eterized particle distributions presented in eqs. (|19p and 
We show explicitly how the parameterizations can 
be used to study correlations between the energy and 
outgoing angle of secondary particles for a general scat- 
tering geometry, viz. two protons with different energies 
colliding at an arbitrary angle. We also present an ex- 
ample in which the energy spectrum of photons due to 
decay of 7r° mesons created in a pp interaction is derived 
from the parameterized ir^ distribution. For clarity, we 
consider only tt'^ mesons in this section, but the presented 
methods are applicable to all pions and kaons. 
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A. Head-on proton — proton collision 

1. Full secondary particle distribution 

First, we consider the energy and rapidity distribution 
of tt'^ mesons created in a collision of an energetic proton 
p with a fixed-target proton q. This is the scattering ge- 
ometry for which the parameterizations presented in this 
work are derived. We denote the Lorentz frame corre- 
sponding to this scattering geometry with K throughout 
this section. In figure (H we present the 7r° distribution 
for incident proton energy Ep = 10^ GeV. We observe 
that the energy and rapidity of the secondary pions are 
strongly correlated: pions with higher energies are emit- 
ted closer to the direction of the incoming energetic pro- 
ton (corresponding to higher values of the rapidity y; see 
eq. dH)). 

Next, we consider two protons p and q colliding head- 
on with energies E'^ and E'^ which defines the reference 
frame K' . Without loss of generality, we take the protons 
to be moving in the z' direction. The secondary particle 
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FIG. 7: Energy spectra of tt" mesons created in a pp collision 
and of the resulting gamma rays. Top panel: lab frame, cor- 
responding to a proton with energy — 10*' GeV colliding 
with a proton at rest; bottom panel: center-of-mass frame, 
corresponding to two protons with energy 730 GeV colliding 
head-on. 
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distribution in this frame is given by 
' e' cosh^ y 



n'(e',y') 



e cosh y' 



n{e, y) , 



(25) 



and the invariance of 
n(e, y) is the particle dis- 



which follows from eq. ((T^ 
en{px,Py,Pz) In this equation, 
tribution in the frame K which is parameterized in eq. 
(fro|l . Note that eq. (pS)) is only valid if the frames K 
and K' are connected with a single Lorentz boost in the 
z (z') direction, i.e. for protons colliding head-on along 
the z' axis in the K' frame. 

As a concrete example, we consider two protons with 
equal energies E'^ — E'^ — 730 GeV. In this case, K' co- 
incides with the center-of-mass (COM) frame for a colli- 
sion between a proton with energy Ep = 10^ GeV and a 
proton at rest. In particular, this means that the center- 
of-mass energy and the secondary multiplicities are 
identical for the scattering geometries in the frames K 
and K'. 

In figure [6l we show the tt^ energy and rapidity distri- 
bution after the collision in the COM frame K' . In this 
frame, the scattering geometry is invariant under the in- 
terchange of the two protons so that the secondary par- 
ticle distribution is symmetric under the transformation 
y —y. It is observed from the figure that this is indeed 
the case for the distribution derived from the parameter- 
ization presented in this work. This is an a posteriori 
verification of the parameterization, which is derived in 
the lab frame without considering this symmetry. 



2. Energy spectrum of secondary particles and decay 
products 

In figure [71 we show the secondary tt" energy spectra 
for the scattering geometries associated with the K and 
K' frames, together with the gamma-ray energy spectra 
resulting from the decay tt'^ — > 77. The decay spectrum 
'T-^(e^) is related to the pion spectrum n{t) as follows 



(see, e.g., Ref. [26 



,(e^) = 2^ 



n(e) 



: de , 



(26) 



where n(e) is the tt" energy spectrum. Because this for- 
mula is valid in all frames, n and e may be replaced by 
n' and e' to derive the gamma-ray energy spectrum from 
the pion energy spectrum in the K' frame. 



B. Proton — proton collision at an arbitrary angle 

In this section, we consider two protons with energies 
i?p and E'^ that collide at an arbitrary angle. Without 
loss of generality, we take proton p to be moving along 
the x' axis in the +x' direction and proton q to be moving 
in the x' - y' plane at an angle 4>'p with respect to the x' 
axis. 




3/2 71 



FIG. 8: Polar plot of the tt" distribution n'(e',0^) as a func- 
tion of the azimuth angle (f>'^ after a collision of a 10* GeV 
proton with a 10^ GeV proton at an angle (f>'q — (3/4)7r. We 
plotted the distribution for pion energies e' = 5 GeV (solid 
line), e' = 1 GeV (dotted line) and e' = 0.5 GeV (dashed 
line). 



We parameterize the distribution of secondary tt^ 
mesons created in this interaction with the pion energy 
e', the zenith angle 9'^ (with respect to the z' axis) and 
the azimuthal angle (f)'^ (in the x' - y' plane). The pion 
momentum is thus expressed as follows; 



k' 



\k'\ sin 9'^ cos (f)'^ ; 



l^'l sine; sin 
I A?' I cos 61; 



(27a) 
(27b) 
(27c) 



where c\k'\ — e'"^ — m'^c'^. In the following, we derive 
the secondary angular distribution in the scattering 
plane and the 7r° energy spectrum. 



1. Secondary angular distribution in the scattering plane 

The pion distribution in the frame K' is derived from 
the parameterization in the fixed-target frame K by 
Lorentz transformations. The frames K' and K are con- 
nected by a Lorentz boost to the rest-frame of proton p, 
followed by a rotation to align the incoming proton q with 
the z axis. The number of secondary pions with energy 
and angles in the range (e' . . . e' + de') x (9'^ . . . 9'^+d9'^) x 
((/.;... cf)'^+d(j)'^) is equal to n'(e', 61;, 0^) sin 9'^d9'^d(f>'^de' , 
where 



n'{e\9'^,cl,'^) 



_ \ n{e,y) 
mld^ -I- e2(l - tanh^ y) ) 27r 



(28) 



In this formula, e and y are the pion energy and rapidity 
in the K frame, respectively, and n(e, y) denotes the pion 



11 




10' 10^ 
e (GeV) 



FIG. 9: Energy spectra of secondary vr" mesons created in 
a collision of a 10* GeV proton with a lO'^ GeV proton for 
three different incident angles 9'^. Also shown is the angle- 
averaged spectrum (see text). For numerical reasons we only 
plot the energy spectrum for head-on collisions at energies 
e > 10^ GeV. We have verified with PYTHIA simulations that 
this part of the spectrum is independent of the incident angle 
between the protons. 




FIG. 10: Polar plot of the tt" distribution n'(e',0^) as a func- 
tion of the zenith angle 9'^ after a collision of a 10* GeV pro- 
ton with an isotropic distribution of 10^ GeV protons. The 
meaning of the lines is the same as in figure [S] 



energy and rapidity distribution which is parameterized 
in eq. (HH). 

In figure [51 we present the distribution of secondary tt'^ 
mesons with respect to the azimuthal angle (/)^, i.e., 

n'(e'>;)= / n'(e',0;>;)sin0;d(?;, (29) 
Jo 

for different values of the pion energy e'. In producing 
this figure, we have chosen incident proton energies E'^ — 
10* GeV and E'^ = 10^ GeV and incident angle ip'^ = 
(3/4)7r. As can be seen from the figure, the pions are 
produced mostly in the direction of the incident protons. 
The degree of collimation is correlated with the energy: 
for pion energies e' below a few GeV, where the pion 
spectrum is highest (see fig. ^ , the pion direction can be 
significantly different from the direction of the colliding 



protons. At energies above a few GeV, the angle of the 
outgoing pion is typically within a few degrees of the 
direction of one of the colliding protons. We have verified 
that this result holds for all secondary pions and kaons. 



2. Secondary energy spectrum 

In figure [3] we present pion energy spectra (integrated 
over pion angles) resulting from a collision of two pro- 
tons with energies E'^ = 10* GeV and E'^ = 10^ GeV 
for different values of the proton collision angle 0^ . For 
comparison, we also show in this figure the pion spectrum 
averaged over incident proton angles (see below). 

While the energy spectrum at high energies is inde- 
pendent of the incident proton angle, there are signifi- 
cant differences at low energies. For small values of the 
incident proton angle 0^, i.e. close to a tail-on collision, 
the low-energy part of the spectrum is suppressed as ex- 
pected. 



C. Isotropic distribution of target protons 

In this section, we consider a single high-energy proton 
p with energy E'^ that interacts with an isotropic distribu- 
tion (in three dimensions) of mono-energetic low-energy 
protons q with energy E'^. We derive the distribution of 
secondary pions with respect to the angle between the 
high-energy proton and the pion, as well as the energy 
spectrum. For an isotropic distribution of target protons, 
the resulting pion distribution does not depend on the 
azimuthal angle around the direction of the high-energy 
proton. In order to keep the former definition of pion an- 
gles (eqs. (|27p ) we consider in this section a high-energy 
incident proton that moves along the z' axis in the +z' 
direction. With this choice, the zenith angle between the 
high-energy proton and the pion is equal to 9'^. 

The momentum of proton q is expressed in terms of 
angles in the same way as the pion momentum in eqs. 
([?7)) : the angle 0'^ denotes the zenith angle with respect 
to the z'-axis and 0^ denotes the azimuth angle with 
respect to the x' axis in the x' - y' plane. 



1. Zenith angle distribution of secondary pions 

The secondary pion distribution, averaged over the in- 
coming angles of low-energy protons q, is given by the 
following expression: 



de'dcos 6: 
1 



d^a' 



<ici de'd cos e'^d(j)'^ 



(30) 



where N is the total number of created pions, cr'^^^^ is 
the inelastic pp cross section and a is the inclusive cross 
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section to detect a particle of a given species assuming 
an ideal detector (cf. section UlI B[) . In this section, we 
use a bar to indicate that a quantity is averaged over the 
incoming angles of low-energy protons q. 

For clarity we assume in this section that both protons 
are very energetic, so that we may take (3'p — f3'^ = 1. The 
averaged inelastic cross section is then equal to 



de'sine'Jl - cos0')ainoi{s{e')) . (31) 



In this equation, cri„ci depends on the proton angle 
through the ccnter-of-mass energy ^/s, where 



2mpC 



(32) 



The dependence of the inelastic cross section on s is ex- 
pressed in eqs. ^ and Q. 

For given values of the proton angles 6'^ and 0'^, the 
differential inclusive cross section and the secondary par- 
ticle distribution are related as follows: 



de'd COS. 6' 



(l-cos0;)a;„,i(s(0;)) 
xn'(e',0;,(/.;;0' (/.'), (33) 



where we have explicitly written the dependence of the 
pion distribution n' on the proton angles 9'^ and (j)'^. The 
total inclusive cross section tr' is obtained by integrat- 
ing eq. (l33|) over the outgoing pion angles and averaging 
over the incident proton angles. The resulting pion dis- 
tribution is homogeneous in the (j)'^ variable. We use this 
rotational invariance to replace the integral over (j)'^ with 
a factor 2tt and choose the value 0^ = to find: 

ct' = de'd9'^d<j)'^ sin 9'^ J d9'^ sin 61^ ( 1 - cos 6*^ ) 

X a[^,M9'^)) n'(e', 0;, (p'^-9'^, 0^ = 0) , (34) 

where the integrals cover the full phase space. The pion 
distribution with respect to the pion energy e' and scat- 
tering angle 9'^ is defined as: 



«'(e',e;) 



d'^N 
de'd9L 



f'^'" d^ N 

^'^We'dcose;^^; 



Using eqs. ([30|) and (p4| . we find that 



n'{e',9'^)^ 



smf 



2 3':. 



d9' sin9'{l~cos9')a[^,,{s{9')) 



inol "'0 
2Tr 



/)>'(£', 0;,(/.;;0'(/.'=O), (36) 



where a[^-^^y is defined in eq. (|3T 



'incl 

In figure [TUl we show the distribution n'{t',9'^) as a 
function of 9'^ for three different values of the pion energy 
e'. In producing this figure, we have considered a collision 
of an energetic proton with energy E'^ = 10^ GeV with 
an isotropic distribution of mono-energetic protons with 



energy E'^ = 10^ GeV. From the figure we observe that 
pions with higher energy are coUimated stronger within 
the direction of the incoming proton, as expected. We 
have verified that this holds for all secondary mesons. In 
the maximally forward direction, i.e., near 9'^ — 0, the 
distribution decreases because the available phase space 
is proportional to sin 9'^ . 



2. Energy spectrum of secondary pions 

The secondary 7r° energy spectrum for the interaction 
of a 10"* GeV proton with the distribution of 10^ GeV 
protons is expressed as 



n'(0= / d9'^n'{e',9'^), 
Jo 



(37) 



where n'(e',6'^) is given in eq. ((36|) . We show in fig- 
ure [9] the energy spectrum averaged over the incoming 
angles of the low-energy protons q. We find that the av- 
eraged spectrum is very close to the spectrum resulting 
from a collision of a 10^ GeV proton with a 10^ GeV pro- 
ton with incident angle 9'^ ~ (5/8)7r, i.e. in the forward 
direction but not head-on. Qualitatively, this is as ex- 
pected because the cross section cr[^^cii^{9'g)) and the flux 
factor (1 — cos 9'^) are largest for head-on collisions while 
the phase-space volume factor sin 9^ suppresses head-on 
collisions. 



VII. EXTRAPOLATION TO THE HIGHEST 
COSMIC-RAY ENERGIES 

The parameterizations presented in section|V]are based 
on simulated pp collisions for incident proton energies 
10^ GeV < Ep < 10^ GeV, where experimental data is 
available to verify the resulting particle distributions and 
multiplicities. Gosmic-ray observations suggest that the 
maximum proton energy that can be generated in astro- 
physical proton accelerators may be as high as 10^^ GeV. 
Thus, in order to account for interactions of the highest 
energy protons, the parameterizations presented in this 
work need to be applied in a region where they cannot 
be directly tested. In this section, we compare the high- 
energy behavior of the parameterizations derived in this 
work with an extrapolation of existing data and with 
theoretical models. Extrapolations of experimental data 
as well as theoretical models for incident proton ener- 
gies Ep > 10^ GeV are available predominantly for the 
charged multiplicity, due to the availability of experimen- 
tal data at lower energies. Therefore, we focus in this 
section on the charged multiplicity contained in the pa- 
rameterizations presented in this work. 

The charged multiplicity is dominated by pions, hence 
we estimate the charged multiplicity from the parameter- 
ized charged pion distributions. We derive the charged 
pion multiplicity A^^± — + M^^- by integrat- 

ing eq. (|19p over energy and rapidity. To account for 
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FIG. 11: Charged multiplicity as a function of incident proton 
energy. The solid line shows an extrapolation of eq. ([Sj; the 
dotted line shows the multiplicity estimated from the parame- 
terized charged pion distributions presented in this work; the 
dashed and dash-dotted lines show the minimum and maxi- 
mum multiphcities given in eqs. (|38p. The data is taken from 

Refs. [33, SIM Ea El- 



charged particle creation due to decay processes and for 
the contribution of other charged particles, we estimate 
the charged multiplicity with M^^' = 2 + lA7Ml±, 
where the numerical value 1.47 is found by comparing 
and Mch at the proton energies considered in our 
simulations. The leading term 2 accounts for the num- 
ber of outgoing protons for low secondary multiplicities 
(corresponding to low center-of-mass energies). 

Using experimental data at low energies, Engel [49| has 
found that the charged multiplicity should increase faster 
than log(s) but not as fast as s^, where 0.1 < p < 0.3, at 
high energies. In order to compare our results with these 
limiting cases, we have re-derived^ the explicit functional 
form based on the two data points with highest energy 
MM: 



7.0+ 1.4 s 



17 logs; 
0.22 



(38a) 
(38b) 



where s is expressed in units of GeV^ . In figure [TT] we 
show the charged multiplicity estimated from the param- 
eterizations presented in this work, together with the 
minimum and maximum values of the multiplicity given 
in eqs. ([35]). Also shown is an extrapolation of the ap- 



proximation M^^ given in eq. ©. We observe that the 
charged multiplicity estimated from our parameteriza- 
tions increases faster than the extrapolation of but 
is well within the limits derived by Engel [4^. We thus 



^ The explicit form of these functions was not given in the original 
work [49| . The numerical value of 0.22 is chosen for comparison 
with fig. 8 of Ref. [H]. 



FIG. 12; Comparison of parameterizations of the 7r° energy 
spectrum en(e) for incident proton energy Ep — 10^ GeV. 
The energy range is chosen for comparison with fig. 1 of Ref. 

a. 



conclude that the high-energy behavior of the parame- 
terizations presented in this work is consistent with the- 
oretical expectations. 

Applying the parameterizations presented in this work 
to very high proton energies Ep > 10^ GeV should be 
done with caution because there is no experimental data 
to directly test the results. Nevertheless, we estimate 
from figure [11] that the uncertainty in the charged parti- 
cle multiplicity at the highest energies Ep ~ 10^^ GeV is 
less than a factor 2. While lacking experimental data on 
the full particle distributions at high energies, we find a 
strong similarity in the particle distributions for different 
proton energies within the three orders of magnitude in 
energy range considered in this work. This suggests that 
the shape of the distributions does not change signifi- 
cantly at very high energies so that the parameterizations 
derived in this work can be applied with some confidence 
to pp collisions with proton energies Ep > 10^ GeV. 



VIII. DISCUSSION 

A. Comparison with previous work 

We have verified that the parameterizations of the pion 
and charged kaon distributions presented in this work are 
similar to those by Badhwar, Stephens and Golden [2^ 
and Stephens and Badhwar 26] at the lowest energies 
considered in this work, ^/s = 45 GeV. The difference 
between our parameterization of the neutral pion distri- 
bution and that of Blattnig et al. [23| is larger. These au- 
thors provide an accurate fit to the particle distribution 
at large transverse momentum. However, the number of 
particles in this region is very small and we find that the 
parameterization does not reproduce the total multiplic- 
ity correctly for center-of-mass energy ^/s = 45 GeV. 
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In figure [121 we present a comparison of the param- 
eterization of the 7r° energy spectrum presented in this 
work with two parameterizations by Kelner, Aharonian 
and Bugayov l24| . These parameterizations are based on 
Monte Carlo results generated with QGSJET and SYBILL 
instead of PYTHIA. We observe from the figure that our 
parameterization is closer to the QGSJET fit at intermedi- 
ate energies and closer to the SYBILL fit at high energies. 
The differences between the energy spectra described by 
the three parameterizations are up to ^^10% for interme- 
diate pion energies. This is larger than the fit inaccu- 
racy (see section [V]) , which suggests that a more precise 
description of the energy spectra and particle distribu- 
tions requires a better theoretical understanding of the 
pp physics, in particular of the fragmentation process, 
rather than more accurate fits. 

The incomplete understanding of the collision physics 
is reflected in the availability of different tuning mod- 
els for PYTHIA. Several tuning models for PYTHIA version 

6.2 can be found in the literature, such as CDF tune A 
(Field^°; see also Ref. [1^) or the tuning proposed in 
Ref. [5l|. In this work, we have used PYTHIA 6.3 with 
default values. A preliminary tuning for PYTHIA version 

6.3 is presented in Ref. ^52i |. It is our intention to adopt 
such a tuning after it has been extensively tested against 
the available data. In this work, we have not studied 
the effect of changing PYTHIA parameters, but we expect 
the differences in the resulting particle distributions to 
be comparable to the those resulting from the use of dif- 
ferent Monte Carlo event generators (see fig. [T^ . 



B. Astrophysical applications 

In section IVTl we demonstrated through explicit exam- 
ples how the parameterizations presented in this work 
can be used to study particle creation in collisions of 
protons with different energies and an arbitrary incident 
angle. We found that the energy of secondary tt'~' mesons 
and the degree of collimation around the direction of the 
colliding protons are strongly correlated (see figures [5] 
and [TU)) . We have verified that this holds for all sec- 
ondary mesons. We expect that this has interesting ob- 
servational consequences for neutrinos and gamma rays 
created in pp collisions in astrophysical sources. 

In the early prompt emission of GRBs, the optical 
depth for pp interactions can be larger than a few. The 
mechanism responsible for the dissipation of the fireball 
energy may accelerate a substantial fraction of the pro- 
tons contained in the fireball to high energies. Subse- 
quent pp collisions give rise to high-energy neutrinos and 
gamma rays. Because the high-energy secondary mesons 
are collimated within the direction of the energetic pro- 



The CDF tunc A is available at the website 
http : //www.phys .uf 1 . edu/'rf ield/cdf /tunes/py_tuneA.html. 



ton, the energy and angular distribution of these secon- 
daries depends on the distribution of the high-energy pro- 
tons. Therefore, the resulting neutrino and gamma-ray 
signals may contain information about the details of the 
acceleration mechanism. We note however that in this 
scenario both interacting protons are moving toward an 
observer with ultra-relativistic velocities so that all par- 
ticle distributions are collimated by relativistic beaming. 
This implies that also low-energy secondaries will be col- 
limated in the observer frame and it will be difficult to 
extract information from the resulting signals. Never- 
theless, the angle-energy correlations found in this work 
may make it possible to extract information about the 
acceleration mechanism. 

The effect of energy-dependent collimation is partic- 
ularly important for observations of secondary particles 
produced in pp collisions where the target proton is non- 
relativistic in the observer frame. A very interesting sce- 
nario in this respect is that of 'failed GRBs' [18, 19] (see 
alsoRefs. [5^[5J|)- It is possible that the mechanisms as- 
sociated with the early phases of a developing GRB may 
be present in a large fraction of supernovae, but only 
lead to an observed GRB under special circumstances. 
For example, it may be the case that the formation of a 
fireball is quite a common phenomenon but that a large 
fraction of fireballs has insufficient energy to traverse the 
pre-burst stellar environment. However, if shocks form 
at a sub-stellar radius, protons can be accelerated and 
collide with target protons (nuclei). The resulting neu- 
trinos are likely the only signals that could indicate the 
existence of such a class of failed (dark) GRBs. In this 
scenario, the energy and flux of the neutrinos that reach 
the earth depend strongly on the collimation of the sec- 
ondary neutrinos created in the pp collisions. 



IX. CONCLUSION 

In this paper we have studied the creation of secondary 
pions and kaons in energetic pp collisions using the event 
generator PYTHIA. We considered an incident proton with 
energy lO'^ GeV < Ep < 10® GeV colhding with a target 
proton at rest. This corresponds to center-of-mass en- 
ergy 43 GeV < < 1-4 X 10^ GeV. The key result is 
the parameterization of the energy spectra (eq. (|16[) ) 
and the energy and rapidity distributions (eqs. (|19p and 
of secondary pions and kaons. Applications of these 
parameterizations are presented in section IVTl In section 
IVIIi we have argued that the results can be applied to 
pp interactions for protons with energies Ep > 10® GeV. 
At the highest CR energies, Ep ~ 10^^ GeV, we have 
estimated the uncertainty in the parameterizations to be 
within a factor 2. 

We have parameterized the particle distributions of 
meta-stable pions and kaons, as opposed to stable decay 
products, because this captures the essential properties 
of the pp interaction without making any assumptions 
about the importance of pion and kaon energy loss prior 
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to decay (for concrete implications of pion decay in an 
astrophysical context, see e.g. Refs. [2^ HO])- Energy 
spectra and full particle distributions of neutrinos and 
gamma rays are derived from the parameterizations pre- 
sented in this work in a straightforward manner. 

The energy and rapidity distributions fully describe 
the kinematics of the secondary mesons, so that the de- 
rived parameterizations contain all correlations between 
energy and angle of the outgoing particles and can be 
applied to a general scattering geometry, viz. two pro- 
tons with different energies that collide under an arbi- 
trary angle. Hence the parameterizations are not limited 
to collisions where one proton is at rest, which opens a 
wealth of astrophysical applications. 

We presented in section I VII several applications of the 
derived parameterizations. We derived the gamma-ray 
spectrum resulting from tt" decay after a pp collision (see 
figure [7]) and we studied correlations between the energy 
and outgoing angle of secondary tt*' mesons produced in 
a collision of two protons with different energies and ar- 
bitrary incident angle (see figure [5]). These results can be 
used for a detailed study of pp interactions in the early 
prompt emission of GRBs and in the interaction of a de- 
veloping GRB with its surroundings (see section IVIIIBp . 
A particularly interesting possibility is the existence of 
a class of developing GRBs for which the fireball has 
insufficient energy to traverse the prc-burst stellar envi- 
ronment. If shocks are formed at a sub-stellar radius, 



these will accelerate protons that collide with target pro- 
tons and create neutrinos. The fluence and energy of 
neutrinos that reach the earth depend sensitively on the 
correlations between the energy and outgoing angle of 
the secondary mesons that are created in these pp inter- 
actions. The parameterizations presented in this work 
can be used to study this scenario in detail. We aim to 
investigate this in the future. 

The parameterizations presented in this work can be 
extended with proton - neutron and proton - photon in- 
teractions, all of which can be studied with PYTHIA. The 
same holds for the energy spectrum and angular distri- 
bution of primary nucleons (the primary nucleon is the 
outgoing nucleon with the highest energy). This allows 
a more precise study of multiple nucleon - nucleon in- 
teractions. These parameterizations are left for future 
work. 
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APPENDIX A: THE LUND STRING MODEL 

The Lund string model [3^ is an iterative model used 
in PYTHIA to describe meson formation after a hard QCD 
process. In the model, quark-antiquark pairs that are 
created in a hard QCD scattering process form 'strings' 
that are connected through a color flux-tube with energy 
per unit length n. This string breaks into a meson and 
a remainder string that will undergo the same process 
(baryons are generated through a generalization of this 
process). At every step in the iteration, a meson is cre- 
ated with a certain energy and rapidity according to a 
probability distribution. 

The mechanism to break the string is the creation of a 
new quark-antiquark pair through quantum-mechanical 
tunneling. The probability to create a qq pair with mass 
m and transverse momentum px is given by 



V = exp 



(-^(mV+p|,c^)), 



(Al) 



which derives from the Schwinger formula [Sa ]. This 
implies that lighter mesons are created more prolifically 
than heavier mesons and that the probability to create 
a meson falls off exponentially with increasing pT- Af- 
ter a meson is created, the probability that it carries a 
fraction z of the string's E + pz is determined by the 
so-called fragmentation function [sil, .32] . Together with 
eq. (|A1[) , this fragmentation function determines the sec- 
ondary particle distributions. Free parameters within the 
model are adjusted to reproduce experimental data. A 
detailed description can be found in Refs. [sills^. 
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TABLE III: Numerical values of the energy spectrum fit parameters pij. 





7r+ 


TT 


7r° 


K+ 


K- 


K" 




Poo 


-0.497 


-0.501 


-0.456 


-1.23 


-1.46 


-1.29 


-1.50 


Poi 


0.0934 


0.0934 


0.0950 


0.0657 


0.101 


0.0690 


0.101 


Pw 


-0.140 


-0.128 


-0.142 


-0.147 


-0.109 


-0.142 


-0.118 


Pll 


0.0131 


0.0118 


0.0135 


0.0161 


0.00865 


0.0154 


0.0101 


P2Q 


-0.455 


-0.437 


-0.457 


-0.00411 


-0.0577 


-0.00717 


-0.0567 


P21 


0.495 


0.494 


0.494 


0.493 


0.491 


0.493 


0.491 


P30 


-2.06 


-0.945 


-1.49 


-0.989 


-1.22 


-1.03 


-1.25 


P40 


-0.896 


-1.03 


-0.981 


-0.345 


-0.164 


-0.294 


-0.169 


P50 


1.11 


0.963 


1.01 


0.777 


1.04 


0.839 


1.05 


Peo 


0.791 


0.598 


0.723 


-0.235 


-0.279 


-0.272 


-0.272 


P70 


37.7 


15.3 


22.1 


42.7 


18.6 


33.8 


21.2 


P80 


7.69 


7.23 


8.53 


12.0 


4.23 


10.1 


4.07 



TABLE IV: Numerical values of the energy and rapidity distribution fit parameters pij (refitted modified energy spectrum) 
and Qij. A hyphen indicates that the parameter is not used in the parameterization. 





7r+ 


TT 


7r° 


K+ 


K- 


K° 




Poo 


-0.474 


-0.461 


-0.420 


-1.03 


-1.20 


-1.13 


-1.37 


Poi 


0.0846 


0.0796 


0.0821 


-0.00299 


0.0377 


0.0168 


0.0501 


pio 


-0.115 


-0.124 


-0.118 


-0.0375 


-0.0291 


-0.0604 


-0.0579 


Pll 


0.0102 


0.0117 


0.0107 


0.00396 


-0.000110 


0.00621 


0.00351 


P20 


-0.560 


-0.604 


-0.598 


-0.835 


-0.606 


-0.655 


-0.467 


P21 


0.497 


0.496 


0.497 


0.494 


0.497 


0.497 


0.497 


P30 


-1.15 


-0.641 


-0.815 


-0.742 


-0.845 


-0.788 


-0.917 


P40 


-1.03 


-1.11 


-1.17 


-0.167 


-0.155 


-0.237 


-0.176 


P50 


1.12 


0.980 


0.987 


0.716 


0.934 


0.840 


1.01 


Pea 


0.962 


0.891 


0.954 


1.08 


0.597 


0.789 


0.371 


P70 
















P80 


6.98 


6.93 


7.45 










goo 


-0.167 


-0.149 


-0.161 


0.539 


0.363 


0.405 


0.382 


qoi 


0.0497 


0.108 


0.0737 


0.222 


0.228 


0.149 


0.195 


qio 


0.668 


0.668 


0.637 


0.889 


1.06 


0.997 


1.12 


?ii 


0.329 


0.328 


0.307 


0.523 


0.673 


0.612 


0.719 


qi2 


0.116 


0.0806 


0.107 


0.227 


0.328 


0.216 


0.298 


gi3 


-0.162 


-0.154 


-0.144 


-0.304 


-0.141 


-0.306 


-0.155 


520 


2.03 


2.05 


2.16 


0.902 


0.695 


0.676 


0.579 


521 


-0.0577 


-0.0654 


-0.0704 


-0.0694 


-0.0648 


-0.0525 


-0.0527 


522 


0.247 


0.233 


0.216 


0.185 


0.226 


0.242 


0.261 


523 


0.665 


0.381 


0.556 














530 


1.04 


1.24 


1.37 


0.319 


0.198 


0.238 


0.184 


531 


3.94 


4.51 


4.92 


1.16 


1.17 


0.951 


1.08 


532 


-1.37 


-1.54 


-1.65 


-0.597 


-0.699 


-0.559 


-0.684 



